source("utils.R")
source("gamm4_utils.R")
load_all_pkgs()
theme_set(theme_bw())
zmargin <- theme(panel.spacing=grid::unit(0,"lines"))
zmarginx <- theme(panel.spacing.x=grid::unit(0,"lines"))
library('pander') ## for tables
sapply(pkgList,function(x) format(packageVersion(x)))
## lme4 gamm4 bbmle broom.mixed brms
## "1.1.21.9002" "0.2.5" "1.0.23.1" "0.2.4.9000" "2.11.1"
## data.table lattice gridExtra ggplot2 viridis
## "1.12.8" "0.20.38" "2.3" "3.2.1" "0.5.1"
## plotly cowplot ggstance plyr dplyr
## "4.9.1" "1.0.0" "0.3.3" "1.8.5" "0.8.4"
## tidyr tibble remef r2glmm raster
## "1.0.2" "2.1.3" "1.0.6.9000" "0.1.2" "3.0.12"
## rgdal fields plotrix sp
## "1.4.8" "10.3" "3.7.7" "1.3.2"
comb_out <- function(p,fn,...) {
print(p)
htmlwidgets::saveWidget(ggplotly(p,...),fn)
}
best_models <- readRDS("bestmodels_gamm4.rds")
ecoreg <- readRDS("ecoreg.rds")
ecoreg <- ecoreg
## names are the names of the variable in the data set;
## the character strings are the names to put on the axis
pred_names <- c(NPP_mean="NPP g*C/m2*year",
NPP_cv_inter="CV of NPP",
Feat_mean="Feat (% of NPP)",
Feat_cv_inter = "CV of Feat")
taxon_names <- paste(c("Birds","Amphibians","Mammals"),
"N sp/km2")
names(taxon_names) <- c("mbirds","mamph","mmamm")
## plot raw values of diversity for a given taxon/predictor
do_plots_1 <- function(taxon="mbirds",pred="NPP_mean",
s_names=taxon_names,
p_names=pred_names) {
ggplot(ecoreg,aes_string(x=pred,y=taxon,colour="biome")) +
geom_point() +
theme(legend.position="none") +
labs(y=s_names[taxon], x=p_names[pred])
}
## generate all combinations
orig_plots <- list()
for (tt in names(taxon_names)) {
orig_plots[[tt]] <- list()
for (pp in names(pred_names)) {
orig_plots[[tt]][[pp]] <- do_plots_1(tt,pp)
}
}
do.call(grid.arrange,c(orig_plots[["mbirds"]],list(nrow=2)))
## names are the names of the variable in the data set;
## the character strings are the names to put on the axis
## this time using log-scaled etc.
sc_pred_names <- c(NPP_log_sc="log NPP g*C/m2*year",
NPP_cv_sc="scaled CV of NPP",
Feat_log_sc="log Feat (% of NPP)",
Feat_cv_sc = "scaled CV of Feat")
sc_taxon_names <- paste("log",taxon_names)
names(sc_taxon_names) <- paste0(names(taxon_names),"_log")
## do all combinations
sc_plots <- list()
for (tt in names(sc_taxon_names)) {
sc_plots[[tt]] <- list()
for (pp in names(sc_pred_names)) {
sc_plots[[tt]][[pp]] <- do_plots_1(tt,pp,
p_names=sc_pred_names,
s_names=sc_taxon_names)
}
}
## for a given taxon, draw unscaled plots in row 1 and scaled plots in row 2
both_plots <- function(taxon) {
do.call(grid.arrange,c(orig_plots[[taxon]],
sc_plots[[paste0(taxon,"_log")]],
list(nrow=2)))
}
FIXME: y-axis names
I generated several “lists” of plots (using plotfun) with the lapply function following Ben’s code. In each case, I changed the xvar and auxvar to plot the 4 top most variables for each taxa (aside from NPP)
taxon <- "mbirds_log"
plotfun(best_models[[taxon]],
respvar=taxon,
xvar='Feat_log_sc',auxvar=NULL,data=ecoreg,backtrans=TRUE,log="xy")
## partial residuals plots
remef_plot <- function(taxon="mbirds_log",predvar="NPP_log",
auxvar=NULL,title=NULL) {
m <- best_models[[taxon]]
if (is.null(title)) {
title <- if (is.null(auxvar)) predvar else {
paste(predvar,auxvar,sep=":")
}
}
pp <- (plotfun(m,xvar=predvar,respvar="partial_res",
auxvar=auxvar,data=ecoreg,
backtrans=TRUE,log="xy")
+ theme(legend.position="none")
+ ggtitle(title)
)
return(pp)
}
## rp1: NPP_log
## rp2: NPP_log:Feat_cv_sc
## rp3: NPP_log:Feat_log
## rp4: Feat_log:NPP_cv_sc
## rp5: Feat_log
## rp6: NPP_cv_sc
## rp7: Feat_cv_sc
rem_predvars <- list("NPP_log_sc",
"NPP_log_sc",
"NPP_log_sc",
"Feat_log_sc",
"Feat_log_sc",
"NPP_cv_sc",
"Feat_cv_sc")
rem_auxvars <- list(NULL,"Feat_cv_sc","Feat_log_sc",
"NPP_cv_sc",NULL,NULL,NULL)
## construct all combinations of partial residuals for each taxon
remef_plots <- list()
for (tt in names(sc_taxon_names)) { ## for each taxon
remef_plots[[tt]] <- list()
for (pp in seq_along(rem_predvars)) { ## for each predictor variable
## cat(tt,pp,"\n")
## cat(tt," ",pp,"\n")
## generate and save the partial residuals plot
remef_plots[[tt]][[pp]] <- remef_plot(tt,predvar=rem_predvars[[pp]],
auxvar=rem_auxvars[[pp]])
## print(remef_plots[[tt]][[pp]])
}
}
ggplot(ecoreg,aes(NPP_mean,mbirds,colour=biome)) + geom_point() + labs(x = "NPP (g*C/m2*year)", y = "Birds (N sp/km2)")
both_plots("mamph")
do.call(grid.arrange,
c(remef_plots[["mamph_log"]][c(2,3,4,5)],
## NPP_log:Feat_csv, NPP_log:Feat_log, Feat_log:NPP_cv_sv, Feat_log
list(ncol=2)))
both_plots("mbirds")
do.call(grid.arrange,
c(remef_plots[["mbirds_log"]][c(5,2,6,7)],
list(ncol=2)))
both_plots("mmamm")
do.call(grid.arrange,
c(remef_plots[["mmamm_log"]][c(7,5,2,4)],
list(ncol=2)))